Statistical analysis using causal inference

Agalychnis lemur

Author

Marcelo Araya-Salas, PhD

Statistical analysis for the paper

In review. Environmental drivers of calling activity in the endangered species Lemur Leaf frog

Purpose

  • Determine the adjustment sets that allow to infer a causal effect of environmental variables on vocal activity
  • Evaluate causal effect of environmental factors on vocal activity of A. lemur with bayesian regression models

Instructions

  • The only input data needed is the file ‘acoustic_and_climatic_data_by_hour.csv’ which has been shared as supplementary data

  • Statistical models can take some time to fit, depending on computing power

  • Many output files are created and re-read in subsequent chunks

Analysis flowchart

This flowchart depicts the data analysis steps described in this report:

flowchart
  A[Define DAG] --> B(24 hours previous rain) 
  A --> C(48 hours previous rain)
  B --> D(Define adjustment sets for each predictor)
  C --> D
  D --> E(Run all models satisfying\nthe back door criterion)
  E --> F(Average posterior probabilities) 
  F --> G(Combine models in a single graph) 

style A fill:#44015466
style B fill:#3E4A894D
style C fill:#3E4A894D
style D fill:#26828E4D
style E fill:#6DCD594D
style F fill:#FDE7254D
style G fill:#31688E4D

Load packages

The code below installs all the necessary packages to run the analyses described in the report:

Code 
if (!require("sketchy", character.only = TRUE))
  install.packages("sketchy")

pkgs <-
  c(
    "remotes",
    "viridis",
    "brms",
    "cowplot",
    "posterior",
    "readxl",
    "HDInterval",
    "kableExtra",
    "knitr",
    "ggdist",
    "ggplot2",
    "lunar",
    "cowplot",
    "maRce10/brmsish",
    "warbleR",
    "ohun",
    "dagitty",
    "ggdag",
    "tidybayes",
    "pbapply"
  )

# install/ load packages
sketchy::load_packages(packages = pkgs)
  
options("digits" = 6, "digits.secs" = 5, knitr.table.format = "html") 

# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE, tidy = TRUE)

# set working directory as project directory or one directory above,
opts_knit$set(root.dir = "..")

Custom functions

Here we create some functions to format data and model outputs:

Code 
adjustment_set_formulas <- function(dag, exposure, required_variable, outcome, effect = "total",
    type = "minimal", formula_parts = c(outcome), latent = NULL, remove = NULL, plot = TRUE,
    ...) {
    if (plot)
        gg <- ggdag_adjustment_set(.tdy_dag = tidy_dagitty(dag), exposure = exposure,
            outcome = outcome, ...) + theme_dag()

    temp_set <- adjustmentSets(x = dag, exposure = exposure, outcome = outcome, effect = effect,
        type = type)

    form_set <- lapply(temp_set, function(x) {
        if (!is.null(remove))
            x <- x[!x %in% remove]
        form <- paste(formula_parts[1], " ~ ", exposure, " + ", paste(x, collapse = " + "),
            if (length(formula_parts) == 2)
                formula_parts[2] else NULL)

        return(form)
    })

    form_set <- form_set[!duplicated(form_set)]

    if (!is.null(latent))
        for (i in latent) form_set <- form_set[!grepl(paste0(" ", i, " "), form_set)]

    # form_set <- sapply(form_set, as.formula)

    names(form_set) <- seq_along(form_set)

    # add formula as attribute
    attributes(form_set)$exposure.formula <- paste(formula_parts[1], " ~ ", exposure,
        if (length(formula_parts) == 2)
            formula_parts[2] else NULL)

    if (plot)
        return(gg) else return(form_set)
}

# Define a function to remove special characters
remove_special_chars <- function(text) {
    # Replace special characters with hyphen
    cleaned_text <- gsub("[^a-zA-Z0-9]+", "-", text)
    # Remove leading and trailing hyphens
    cleaned_text <- gsub("^-+|-+$", "", cleaned_text)
    return(cleaned_text)
}

pa <- function(...) brms::posterior_average(...)

# to get average posterior values from models with different formulas
averaged_model <- function(formulas, data, model_call, ndraws = 1000, save.files = TRUE,
    path = ".", suffix = NULL, cores = 1, name = NULL) {
    if (dir.exists(file.path(path, name))) {
        cat("Directory already existed. Attempting to fit missing models\n")
        cat("Fitting models (step 1 out of 2) ...")
    } else dir.create(path = file.path(path, name))

    cat("Fitting models (step 1 out of 2) ...")
    fit_list <- pblapply_brmsish_int(X = formulas, cl = cores, function(y) {

        # make file name without special characters
        mod_name <- paste0(remove_special_chars(as.character(y)), ".RDS")

        if (save.files & !file.exists(file.path(path, mod_name))) {
            cat("Fitting", y, "\n")
            mc <- gsub(pattern = "formula", replacement = as.character(y), x = model_call)

            mc <- parse(text = mc)

            fit <- eval(mc)

            if (save.files)
                saveRDS(fit, file = file.path(path, mod_name))

        } else {
            cat("Reading", y, "(already existed)\n")
            fit <- readRDS(file.path(path, mod_name))
        }
        return(fit)
    })

    if (length(formulas) > 1) {
        cat("Averaging models (step 2 out of 2) ...")
        average_call <- parse(text = paste("pa(", paste(paste0("fit_list[[", seq_along(fit_list),
            "]]"), collapse = ", "), ", ndraws = ", ndraws, ")"))

        # Evaluate the expression to create the call object
        average_eval <- eval(average_call)

        # add formula as attribute
        attributes(average_eval)$averaged_fit_formulas <- formulas

        rds_name <- if (is.null(suffix))
            file.path(path, paste0(name, ".RDS")) else file.path(path, paste0(suffix, "_", name, ".RDS"))

        if (save.files)
            saveRDS(average_eval, file = rds_name)

        # return draws from average models
        return(average_eval)
    } else cat("No model averaging conducted as a single formula was supplied")
}

# function to convert effect sizes to percentage of change
to_change_percentage <- function(x) {

    (exp(x) - 1) * 100

}

# prints the results of the averaged models
draw_extended_summary <- function(draws, name = NULL, highlight = TRUE, fill = "#6DCD59FF",
    remove.intercepts = FALSE, by = NULL, gsub.pattern = NULL, gsub.replacement = NULL,
    xlab = "Effect size", ylab = "Parameter") {

    # create objects just to avoid errors with ggplot functions when checking
    # package
    position_dodge <- org.variable <- value <- significance <- `l-95% CI` <- `u-95% CI` <- theme <- unit <- NULL

    if (!is.null(by)) {
        levels <- draws[, by]

        unique_levels <- unique(levels)
    } else unique_levels <- ".A"

    # run loop over each level (or just the single 'level')
    results_list <- lapply(unique_levels, function(x) {
        if (identical(unique_levels, ".A")) {
            # keep only betas
            draws <- draws[, grep("^b_", names(draws), value = TRUE)]
        } else {
            draws <- draws[draws[, by, drop = TRUE] == x, grep("^b_", names(draws),
                value = TRUE)]
        }

        # remove intercept betas
        if (remove.intercepts)
            draws <- draws[, grep("^b_Intercept", names(draws), value = TRUE, invert = TRUE),
                drop = FALSE]

        # compute model-averaged posteriors of overlapping parameters
        coef_table <- posterior::summarise_draws(draws, median, ~quantile(.x, probs = c(0.025,
            0.975)), posterior::default_convergence_measures())

        names(coef_table)[3:4] <- c("l-95% CI", "u-95% CI")

        coef_table$value <- coef_table$median
        coef_table$significance <- ifelse(coef_table$`l-95% CI` * coef_table$`u-95% CI` >
            0, "sig", "non-sig")
        coef_table$significance <- factor(coef_table$significance, levels = c("non-sig",
            "sig"))

        if (ncol(draws) > 1)
            sdraws <- stack(draws[, grep("^b_", names(draws), value = TRUE)], ind = "levels") else sdraws <- data.frame(value = draws[, 1], variable = names(draws))

        names(sdraws) <- c("value", "variable")

        sdraws$significance <- sapply(sdraws$variable, function(x) coef_table$significance[as.character(coef_table$variable) ==
            x][1])

        # add level if (!identical(unique_levels, '.A')) { to coef_table
        coef_table$.new.column <- x
        names(coef_table)[ncol(coef_table)] <- "levels"

        # to sdraws
        sdraws$.new.column <- x
        names(sdraws)[ncol(sdraws)] <- "levels"
        # }

        output <- list(coef_table = coef_table, sdraws = sdraws)

        return(output)
    })

    # put both results into a single data frame
    coef_table <- do.call(rbind, lapply(results_list, "[[", 1))
    sdraws <- do.call(rbind, lapply(results_list, "[[", 2))

    if (!is.null(gsub.pattern) & !is.null(gsub.replacement)) {
        if (length(gsub.pattern) != length(gsub.replacement))
            stop2("'gsub.replacement' and 'gsub.pattern' must have the same length")

        for (i in 1:length(gsub.pattern)) {
            sdraws$variable <- gsub(pattern = gsub.pattern[i], replacement = gsub.replacement[i],
                sdraws$variable)
            coef_table$variable <- gsub(pattern = gsub.pattern[i], replacement = gsub.replacement[i],
                coef_table$variable)
        }
    }

    # duplicate variable column
    coef_table <- as.data.frame(coef_table)
    coef_table$org.variable <- coef_table$variable
    sdraws$org.variable <- sdraws$variable

    if (!identical(unique_levels, ".A")) {
        coef_table$variable <- paste(coef_table$levels, coef_table$variable, sep = "-")
        coef_table <- coef_table[order(coef_table$org.variable), ]

        sdraws$variable <- paste(sdraws$levels, sdraws$variable, sep = "-")
    }

    rownames(coef_table) <- coef_table$variable

    fill_df <- data.frame(level = unique_levels, fill = fill)

    if (!identical(unique_levels, ".A"))
        coef_table$fill_values <- sapply(coef_table$levels, function(x) fill_df$fill[fill_df$level ==
            x]) else coef_table$fill_values <- fill[1]

    if (highlight) {
        coef_table$col_pointrange <- ifelse(coef_table$significance == "non-sig",
            "gray", "black")

        # coef_table$fill_values <- ifelse(coef_table$significance == 'no-sig',
        # grDevices::adjustcolor(col = coef_table$fill_values, alpha.f = 0.5),
        # coef_table$fill_values)
        sdraws$significance <- sapply(sdraws$variable, function(x) coef_table$significance[coef_table$variable ==
            x])
    } else {
        coef_table$col_pointrange <- rep("black", nrow(coef_table))
        sdraws$significance <- 1
    }

    pd <- position_dodge(width = 0.05)

    # creat plots
    gg_distributions <- ggplot2::ggplot(data = sdraws, ggplot2::aes(y = org.variable,
        x = value, fill = levels, alpha = if (highlight)
            significance else NULL)) + ggplot2::geom_vline(xintercept = 0, col = "black", lty = 2) +
        ggdist::stat_halfeye(ggplot2::aes(x = value), .width = c(0.95), normalize = "panels",
            color = "transparent", position = pd) + ggplot2::scale_alpha_manual(values = c(0.4,
        0.8), guide = "none") + ggplot2::scale_fill_manual(values = if (!identical(unique_levels,
        ".A"))
        as.vector(coef_table$fill_values) else fill[1]) + ggplot2::geom_point(data = coef_table, position = pd) + ggplot2::geom_errorbar(data = coef_table,
        ggplot2::aes(xmin = `l-95% CI`, xmax = `u-95% CI`), width = 0, position = pd) +
        ggplot2::facet_wrap(~org.variable, scales = "free_y", ncol = 1) + ggplot2::theme_classic() +
        ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"), plot.margin = ggplot2::margin(0,
            0, 0, 0, "pt"), legend.position = "none", strip.background = ggplot2::element_blank(),
            strip.text = ggplot2::element_blank()) + ggplot2::labs(x = xlab, y = ylab,
        fill = "Effect") + theme(panel.spacing.y = unit(0, "null")  # no vertical space between panels
)

    coef_table$variable <- coef_table$significance <- coef_table$value <- NULL
    names(coef_table) <- c("Estimate", "l-95% CI", "u-95% CI", "Rhat", "Bulk_ESS",
        "Tail_ESS")

    coef_table <- coef_table[, c("Estimate", "l-95% CI", "u-95% CI", "Rhat", "Bulk_ESS",
        "Tail_ESS")]

    if (!is.null(name))
        cat("\n\n## ", name, "\n\n")

    html_coef_table <- brmsish:::html_format_coef_table(coef_table, fill = fill,
        highlight = highlight)

    print(html_coef_table)

    print(gg_distributions)
}

1 Read data

1.1 Set working directory

This is the path were the supplementary data is saved and were the output files will be saved. It also creates a directory there in which to save single and averaged models:

Code 
# Set working directory
path <- "PUT_WORKING_DIRECTORY_PATH_HERE"

# create directory for single models
dir.create(file.path(path, "single_models"))

# create directory for average models
dir.create(file.path(path, "averaged_models"))

The file ‘acoustic_and_climatic_data_by_hour.csv’ contains all the data required to run statistical models down below:

Code 
call_rate_hour <- read.csv(file.path(path, "acoustic_and_climatic_data_by_hour.csv"))

Describe environmental variables:

Code 
agg <- aggregate(cbind(temp, prev_temp, HR, rain, rain_24, rain_48, moonlight) ~
    1, call_rate_hour, function(x) round(c(mean(x), sd(x), min(x), max(x)), 3))

agg <- as.data.frame(matrix(unlist(agg), ncol = 4, byrow = TRUE, dimnames = list(c("Temperature",
    "Previous temperature", "Relative humidity", "Night rain", "Rain 24 hours", "Rain 48 hours",
    "Moonlight"), c("mean", "sd", "min", "max"))))
# print table as kable
kb <- kable(agg, row.names = TRUE, digits = 3)

kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))

print(kb)
mean sd min max
Temperature 24.229 2.126 19.472 31.763
Previous temperature 24.097 0.988 20.005 27.111
Relative humidity 90.056 8.464 55.158 99.940
Night rain 0.079 0.457 0.000 10.417
Rain 24 hours 1.441 3.161 0.000 25.620
Rain 48 hours 1.453 3.170 0.000 25.620
Moonlight 0.495 0.355 0.000 1.000

2 Directed acyclical graphs (DAGs)

The code creates the DAGs used for determining the adjustments sets needed to evaluate causality:

Code 
coords <- list(x = c(sc_rain = -0.4, evotranspiration = 0.5, sc_prev_rain = 0.7,
    sc_temp = -0.8, sc_HR = 0, n_call = 0, sc_moonlight = 0.3, hour = -0.5), y = c(sc_rain = 0.4,
    evotranspiration = 0.3, sc_prev_rain = -0.5, sc_temp = 0, climate = 1, sc_HR = -0.6,
    n_call = 0, sc_moonlight = 1, hour = 0.9))

# sc_temp + sc_HR + sc_moonlight + sc_rain + sc_rain_24 + ar(p = 2, time =
# hour_diff, gr = hour sc_temp = temp y meanT = prev_temp

dag_l <- dagify(sc_rain ~ evotranspiration, sc_prev_rain ~ evotranspiration, sc_temp ~
    climate, sc_temp ~ sc_rain, sc_HR ~ sc_rain, n_call ~ sc_HR, n_call ~ hour, n_call ~
    sc_moonlight, sc_moonlight ~ hour, sc_temp ~ hour, sc_HR ~ sc_temp, sc_HR ~ sc_prev_rain,
    sc_HR ~ sc_rain, n_call ~ sc_temp, n_call ~ sc_prev_rain, n_call ~ sc_rain, labels = c(n_call = "Calling\nactivity",
        sc_HR = "Relative\nhumidity", sc_rain = "Current\nrain", sc_prev_rain = "Prior\nrain",
        sc_moonlight = "Moonlight", hour = "Earth\nrotation", sc_temp = "Tempera-\nture",
        evotranspiration = "Evotrans-\npiration", climate = "Climate", latent = c("evotranspiration",
            "climate"), outcome = "n_call"), coords = coords)

tidy_dag <- tidy_dagitty(dag_l)
tidy_dag$data$type <- ifelse(is.na(tidy_dag$data$to), "outcome", "predictor")
tidy_dag$data$type[tidy_dag$data$name %in% c("evotranspiration", "climate")] <- "latent"


dat <- tidy_dag$data
shorten_distance <- c(0.07, 0.07)
dat$slope <- (dat$yend - dat$y)/(dat$xend - dat$x)
distance <- sqrt((dat$xend - dat$x)^2 + (dat$yend - dat$y)^2)
proportion <- shorten_distance[1]/distance
dat$xend <- (1 - proportion/2) * dat$xend + (proportion/2 * dat$x)
dat$yend <- (1 - proportion/2) * dat$yend + (proportion/2 * dat$y)
proportion <- shorten_distance[2]/distance
dat$xstart <- (1 - proportion/2) * (dat$x - dat$xend) + dat$xend
dat$ystart <- (1 - proportion/2) * (dat$y - dat$yend) + dat$yend

tidy_dag$data <- dat


basic_dag <- ggplot(tidy_dag, aes(x = x, y = y, xend = xend, yend = yend)) + scale_color_viridis_d(begin = 0.2,
    end = 0.8, alpha = 0.5) + geom_dag_text(color = "black", aes(label = label, color = label)) +
    labs(color = "Type") + theme_dag() + theme(legend.position = "bottom") + guides(colour = guide_legend(override.aes = list(size = 10))) +
    geom_dag_point(aes(color = type), size = 30) + expand_limits(y = c(-0.67, 1.1)) +
    geom_dag_edges_fan(edge_color = viridis(10, alpha = 0.4)[2], arrow = grid::arrow(length = grid::unit(10,
        "pt"), type = "closed"), aes(x = xstart, y = ystart, xend = xend, yend = yend))

basic_dag

Code 
dag_24 <- dagify(sc_rain ~ evotranspiration, sc_rain_24 ~ evotranspiration, sc_temp ~
    climate, sc_temp ~ sc_rain, sc_HR ~ sc_rain, n_call ~ sc_HR, n_call ~ hour, n_call ~
    sc_moonlight, sc_moonlight ~ hour, sc_temp ~ hour, sc_HR ~ sc_temp, sc_HR ~ sc_rain_24,
    sc_HR ~ sc_rain, n_call ~ sc_temp, n_call ~ sc_rain_24, n_call ~ sc_rain, labels = c(n_call = "Call rate",
        sc_HR = "Relative humidity", sc_rain = "Night Rain", sc_rain_24 = "Previous Rain",
        sc_moonlight = "Moonlight", hour = "Earth rotation", sc_temp = "Temperature",
        evotranspiration = "Evotranspiration", climate = "Climate", latent = c("evotranspiration",
            "climate"), outcome = "n_call"))

dag_48 <- dagify(sc_rain ~ evotranspiration, sc_rain_48 ~ evotranspiration, sc_temp ~
    climate, sc_temp ~ sc_rain, sc_HR ~ sc_rain, n_call ~ sc_HR, n_call ~ hour, n_call ~
    sc_moonlight, sc_moonlight ~ hour, sc_temp ~ hour, sc_HR ~ sc_temp, sc_HR ~ sc_rain_48,
    sc_HR ~ sc_rain, n_call ~ sc_temp, n_call ~ sc_rain_48, n_call ~ sc_rain, labels = c(n_call = "Call rate",
        sc_HR = "Relative humidity", sc_rain = "Night Rain", sc_rain_48 = "Previous Rain",
        sc_moonlight = "Moonlight", hour = "Earth rotation", sc_temp = "Temperature",
        evotranspiration = "Evotranspiration", climate = "Climate", latent = c("evotranspiration",
            "climate"), outcome = "n_call"))

3 Bayesian regression models

3.1 Scale variables and set model parameters

Prepare data for regression models:

Code 
# make hour a factor
call_rate_hour$hour <- factor(call_rate_hour$hour)

# scale and mean-center
call_rate_hour$sc_temp <- scale(call_rate_hour$temp)
call_rate_hour$sc_HR <- scale(call_rate_hour$HR)
call_rate_hour$sc_rain <- scale(call_rate_hour$rain)
call_rate_hour$sc_rain_24 <- scale(call_rate_hour$rain_24)
call_rate_hour$sc_rain_48 <- scale(call_rate_hour$rain_48)
call_rate_hour$sc_moonlight <- scale(call_rate_hour$moonlight)

priors <- c(prior(normal(0, 4), class = "b"))
chains <- 4
iter <- 10000

3.2 Fit all models

This code fits all the models in the adjustment sets across all predictors:

Code 
param_grid <- expand.grid(dag = c("dag_24", "dag_48"), exposure = c("sc_temp", "sc_HR",
    "sc_moonlight", "sc_rain", "sc_rain_24", "sc_rain_48"), effect = c("total", "direct"),
    stringsAsFactors = FALSE)

param_grid$name <- apply(param_grid, 1, paste, collapse = "-")

# remove wrong dags for previous rain
param_grid <- param_grid[!(param_grid$dag == "dag_24" & param_grid$exposure == "sc_rain_48") &
    !(param_grid$dag == "dag_48" & param_grid$exposure == "sc_rain_24"), ]


adjustment_sets_list <- pblapply(seq_len(nrow(param_grid)), cl = 1, function(x) {

    forms <- adjustment_set_formulas(dag = if (param_grid$dag[x] == "dag_24")
        dag_24 else dag_48, type = if (param_grid$effect[x] == "total")
        "all" else "minimal", exposure = param_grid$exposure[x], outcome = "n_call", effect = param_grid$effect[x],
        required_variable = "hour", formula_parts = c("n_call | resp_rate(rec_time)",
            "+ ar(p = 2, time = hour_diff, gr = hour)"), latent = c("evotranspiration",
            "climate"), remove = "hour", plot = FALSE)

    return(forms)
})



names(adjustment_sets_list) <- param_grid$name

param_grid$model <- sapply(seq_len(nrow(param_grid)), function(x) {
    if (param_grid$effect[x] == "direct")
        adjustment_sets_list[[which(names(adjustment_sets_list) == param_grid$name[x])]] else NA
})

param_grid$model <- unlist(param_grid$model)
param_grid <- as.data.frame(param_grid)


param_grid$model[!is.na(param_grid$model)] <- remove_special_chars(param_grid$model[!is.na(param_grid$model)])
param_grid$model <- c("total_effect_temperature_with_rain_24", "total_effect_temperature_with_rain_48",
    "total_effect_humidity_with_rain_24", "total_effect_humidity_with_rain_48", "total_effect_moon_with_rain_24",
    "total_effect_moon_with_rain_48", "total_effect_rain_with_rain_24", "total_effect_rain_with_rain_48",
    "total_effect_previous_rain_24", "total_effect_previous_rain_48", param_grid$model[!is.na(param_grid$model)])

param_grid$exposure.name <- param_grid$exposure
param_grid$exposure.name[grep("temp", param_grid$exposure.name)] <- "Temperature"
param_grid$exposure.name[grep("HR", param_grid$exposure.name)] <- "Relative humidity"
param_grid$exposure.name[grep("moon", param_grid$exposure.name)] <- "Moonlight"
param_grid$exposure.name[grep("rain$", param_grid$exposure.name)] <- "Current rain"
param_grid$exposure.name[grep("rain_24", param_grid$exposure.name)] <- "Previous rain (24h)"
param_grid$exposure.name[grep("rain_48", param_grid$exposure.name)] <- "Previous rain (48h)"

table(param_grid$exposure.name)

write.csv(x = param_grid, file = file.path(path, "direct_and_total_effect_model_data_frame.csv"),
    row.names = FALSE)

3.3 Combined models need to infer causality

The code below takes all models representing an adjustment set for each predictor and average them into a single model fit:

Code 
direct_adjustment_sets_list <- adjustment_sets_list[grep("direct", names(adjustment_sets_list))]

for(i in seq_along(direct_adjustment_sets_list))
pa_comb_mod <-
averaged_model(
    formulas = direct_adjustment_sets_list[[i]],
    data = call_rate_hour,
    suffix = "direct",
    model_call = "brm(formula, data = call_rate_hour, iter = iter, chains = chains, cores = chains, family = negbinomial(), prior = priors, backend = 'cmdstanr')",
    save.files = TRUE,
    path = file.path(path, "averaged_models"), 
    # name = "temperature_with_rain_24",
    cores = 1 
  )


model_call = "brm(formula, data = call_rate_hour, iter = iter, chains = chains, cores = chains, family = negbinomial(), prior = priors, backend = 'cmdstanr')"
formulas <- unlist(direct_adjustment_sets_list)
mod_path <- file.path(path, "single_models")

fit_list <- pblapply_brmsish_int(X = formulas, cl = 1, function(y) {
        # make file name without special characters
        mod_name <-
          paste0(remove_special_chars(as.character(y)), ".RDS")
        
        if (!file.exists(file.path(mod_path, mod_name))) {
          cat("Fitting", y, "\n")
          mc <-
            gsub(pattern = "formula",
                 replacement = as.character(y),
                 x = model_call)
          
          mc <- parse(text = mc)
          
          fit <- eval(mc)
          
          if (save.files)
            saveRDS(fit, file = file.path(mod_path, mod_name))
          
        } 
      })

3.4 Results

Code 
param_grid <- read.csv(file = file.path(path, "direct_and_total_effect_model_data_frame.csv"))

param_grid$files <- file.path(file.path(path, "averaged_models", paste0(param_grid$model,
    ".RDS")))


for (i in unique(param_grid$exposure.name)) {
    Y <- param_grid[param_grid$exposure.name == i, ]

    cat(paste("\n###", i), "\n")
    cat("\n#### Direct effects\n")
    for (e in which(Y$effect == "direct")) {
        if (grepl("24", Y$model[e]))
            cat("\n##### 24 hour previous rain model:\n") else cat("\n##### 48 hour previous rain model:\n")
        extended_summary(read.file = Y$files[e], highlight = TRUE, remove.intercepts = TRUE,
            print.name = FALSE)
        cat("\n")
    }

    cat("\n#### Total effect\n")
    for (w in which(Y$effect == "total")) {
        if (grepl("24", Y$files[w]))
            cat("\n##### 24 hour previous rain model:\n") else cat("\n##### 48 hour previous rain model:\n")

        draws <- readRDS(Y$files[w])

        draw_extended_summary(draws, highlight = TRUE, remove.intercepts = TRUE)
        cat("\n")
        cat("\n###### Summary of single models:\n")

        # print summary
        print(readRDS(gsub("\\.RDS", "_fit_summary.RDS", Y$files[w])))
    }
    cat("\n")
}

3.4.1 Temperature

3.4.1.1 Direct effects

3.4.1.1.1 24 hour previous rain model:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 13808 14519.5 1687235322
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_temp 0.577 0.498 0.656 1 13808.0 14685.6
b_sc_HR 0.314 0.242 0.387 1 18992.0 14519.5
b_sc_rain 0.037 0.002 0.074 1 23964.4 15266.3
b_sc_rain_24 0.096 0.057 0.135 1 23319.6 15914.6

3.4.1.1.2 48 hour previous rain model:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 9460.49 12432.4 293904519
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_temp 0.565 0.487 0.643 1 9460.49 12432.4
b_sc_HR 0.307 0.235 0.379 1 11994.69 13312.6
b_sc_rain 0.029 -0.007 0.066 1 17566.34 16050.0
b_sc_rain_48 0.006 -0.032 0.044 1 17289.44 16222.0

3.4.1.2 Total effect

3.4.1.2.1 24 hour previous rain model:
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_temp 0.308 0.258 0.361 1.043 44.190 893.839
b_sc_rain 0.041 0.004 0.078 1.026 152.336 830.415

3.4.1.2.1.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2542.43 5383.76 629626678
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2238.34 5028.64 2096770826
3 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2320.90 4619.08 1430954384
4 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2292.54 4473.95 1881604582
3.4.1.2.2 48 hour previous rain model:
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_temp 0.305 0.258 0.353 1.005 961.292 983.283
b_sc_rain 0.036 -0.002 0.076 1 1113.522 969.564

3.4.1.2.2.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2529.79 5016.37 1471048663
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2588.68 5844.63 412112948
3 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2622.63 5641.91 41093237
4 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2334.61 4741.38 366114719

3.4.2 Relative humidity

3.4.2.1 Direct effects

3.4.2.1.1 24 hour previous rain model:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 11039.8 13170.1 1554924087
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_HR 0.314 0.240 0.388 1 13006.8 13679.9
b_sc_rain 0.037 0.002 0.075 1 18235.1 15400.0
b_sc_rain_24 0.096 0.057 0.135 1 16966.8 15476.0
b_sc_temp 0.578 0.499 0.657 1 11039.8 13170.1

3.4.2.1.2 48 hour previous rain model:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 8448.99 11500.1 1931567299
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_HR 0.308 0.234 0.380 1 10649.49 12433.6
b_sc_rain 0.029 -0.008 0.067 1 16664.06 15531.5
b_sc_rain_48 0.005 -0.032 0.043 1 15745.61 14835.4
b_sc_temp 0.566 0.487 0.645 1 8448.99 11500.1

3.4.2.2 Total effect

3.4.2.2.1 24 hour previous rain model:
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_HR 0.309 0.240 0.383 1.006 942.064 903.519
b_sc_rain 0.038 0.002 0.076 1.003 995.090 912.254
b_sc_rain_24 0.094 0.057 0.133 1.003 993.441 772.255
b_sc_temp 0.575 0.496 0.655 1.004 1022.397 861.998

3.4.2.2.1.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_HR + sc_moonlight + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2556.01 5127.35 1982612331
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2182.06 4328.43 1554924087
3.4.2.2.2 48 hour previous rain model:
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_HR 0.303 0.230 0.371 1.006 801.527 931.946
b_sc_rain 0.030 -0.008 0.066 1.002 812.238 1071.444
b_sc_rain_48 0.003 -0.034 0.044 1.005 889.913 984.434
b_sc_temp 0.564 0.480 0.644 1.003 821.550 773.876

3.4.2.2.2.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_HR + sc_moonlight + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2326.55 4696.55 846435536
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2166.28 4627.63 1931567299

3.4.3 Moonlight

3.4.3.1 Direct effects

3.4.3.1.1 48 hour previous rain model:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 30283.5 15929.8 1624275037
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_moonlight -0.129 -0.193 -0.064 1 30283.5 15929.8

3.4.3.1.2 48 hour previous rain model:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 30283.5 15929.8 1624275037
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_moonlight -0.129 -0.193 -0.064 1 30283.5 15929.8

3.4.3.2 Total effect

3.4.3.2.1 24 hour previous rain model:
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_moonlight -0.101 -0.17 -0.039 1.04 40.665 750.448

3.4.3.2.1.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2463.32 4856.49 1624275037
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2613.63 5435.51 1331671454
3 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2426.99 4881.26 1714115616
4 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2184.13 4557.11 946500998
5 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2386.65 5083.35 2088023938
6 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2571.23 4971.80 1587538237
7 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2319.91 4626.15 1236984196
8 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2197.61 4746.52 1987527846
9 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2353.89 5263.44 924618251
10 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2409.47 4652.12 1063522846
11 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2407.64 4851.64 259660598
12 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2473.17 4731.37 1625659188
13 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2596.10 5292.48 1259975215
14 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 1998.35 4047.25 192228769
15 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2365.76 4743.96 559459570
16 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2432.09 4364.34 1786523665
3.4.3.2.2 48 hour previous rain model:
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_moonlight -0.105 -0.166 -0.038 1.02 111.506 659.054

3.4.3.2.2.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2534.83 5170.40 169376128
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2664.91 5585.82 731739252
3 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2449.93 4958.79 230289957
4 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2283.15 4266.04 1735286727
5 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2486.75 4561.76 1539258099
6 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2525.61 5053.20 1234320719
7 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2481.99 4862.29 183236870
8 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2296.24 4114.83 239469808
9 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2298.87 4495.74 1126617762
10 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2442.38 4594.40 684931243
11 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 1 (5e-05%) 0 2416.00 4930.19 1093486403
12 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2146.30 4395.84 1670597608
13 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2450.19 4442.71 627408872
14 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2513.43 4849.68 1827601939
15 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2101.51 4415.14 104020182
16 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2495.11 5464.93 598459253
17 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2120.85 4487.97 2041361978
18 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2336.51 4816.06 1759350038
19 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 1 (5e-05%) 0 2453.12 5040.07 927628753
20 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2360.45 5250.70 2028572807
21 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2642.43 5077.10 293339431
22 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2342.88 4887.79 1868346384
23 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2515.33 5176.64 552665161
24 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2253.53 4141.73 1764642823
25 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2634.31 5119.83 505424120
26 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2427.59 5073.57 1761124445
27 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2609.57 5134.76 1073467634
28 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2461.43 4814.97 1419285187
29 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2596.36 5583.40 1232772932
30 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2271.50 5053.09 278133969
31 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 1 (5e-05%) 0 2507.89 4315.91 866725413
32 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_moonlight + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2126.46 4549.18 1885323046

3.4.4 Current rain

3.4.4.1 Direct effects

3.4.4.1.1 24 hour previous rain model:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain + sc_HR + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 10228.7 12456.7 1281527469
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain 0.037 0.001 0.075 1 18004.9 15909.5
b_sc_HR 0.313 0.240 0.386 1 13025.4 13472.6
b_sc_rain_24 0.096 0.059 0.134 1 17444.9 15779.0
b_sc_temp 0.577 0.499 0.655 1 10228.7 12456.7

3.4.4.1.2 48 hour previous rain model:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain + sc_HR + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 7371.36 11882.6 1129340659
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain 0.029 -0.007 0.066 1 14353.91 14479.2
b_sc_HR 0.307 0.235 0.380 1 8890.14 12279.5
b_sc_rain_48 0.005 -0.032 0.043 1 13843.68 14097.9
b_sc_temp 0.566 0.487 0.644 1 7371.36 11882.6

3.4.4.2 Total effect

3.4.4.2.1 24 hour previous rain model:
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain 0.000 -0.034 0.035 1.001 942.79 942.018
b_sc_rain_24 0.076 0.039 0.113 1 1056.52 982.799

3.4.4.2.1.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain + sc_moonlight + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2531.75 4835.39 1596946440
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2369.16 4574.59 1687235322
3.4.4.2.2 48 hour previous rain model:
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain -0.005 -0.042 0.030 1 1024.455 819.064
b_sc_rain_48 0.011 -0.027 0.049 1.009 932.089 975.233

3.4.4.2.2.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain + sc_moonlight + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2647.66 5129.91 1485832077
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2675.26 5099.07 167092432

3.4.5 Previous rain (24h)

3.4.5.1 Direct effects

3.4.5.1.1 24 hour previous rain model:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 11374.2 13445.2 1379661434
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain_24 0.096 0.058 0.135 1 20345.7 15663.3
b_sc_HR 0.314 0.240 0.386 1 14938.4 13445.2
b_sc_rain 0.038 0.001 0.074 1 21502.1 16005.3
b_sc_temp 0.578 0.498 0.656 1 11374.2 13493.9

3.4.5.2 Total effect

3.4.5.2.1 24 hour previous rain model:
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain_24 0.088 0.051 0.127 1.01 883.631 733.445
b_sc_rain 0.043 0.005 0.081 1.001 908.203 446.195

3.4.5.2.1.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2692.30 5204.02 423216939
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2556.12 4635.90 610868624
3 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 1 (5e-05%) 0 2559.09 4449.50 1988187279
4 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2077.52 4134.53 640249946

3.4.6 Previous rain (48h)

3.4.6.1 Direct effects

3.4.6.1.1 48 hour previous rain model:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) n_call | resp_rate(rec_time) ~ sc_rain_48 + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 11112 13056.2 1238572765
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain_48 0.006 -0.031 0.043 1 17211.0 15210.0
b_sc_HR 0.308 0.236 0.380 1 14154.7 13956.7
b_sc_rain 0.029 -0.007 0.066 1 19170.2 15964.2
b_sc_temp 0.566 0.486 0.645 1 11112.0 13056.2

3.4.6.2 Total effect

3.4.6.2.1 48 hour previous rain model:
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sc_rain_48 0.005 -0.033 0.043 1.014 711.511 873.701
b_sc_rain 0.036 0.001 0.073 1 857.372 792.928

3.4.6.2.1.1 Summary of single models:
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2447.89 5823.08 1081516016
2 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 2 (1e-04%) 0 2413.93 4851.86 217409394
3 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2506.36 5173.40 724206553
4 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2641.41 5493.32 509251240
5 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_48 + sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2386.41 4640.41 518390312
6 b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) n_call | resp_rate(rec_time) ~ sc_rain_48 + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) 10000 4 1 5000 0 (0%) 0 2575.04 5332.87 76427694

3.5 Combined results with causal inference estimates

Takes the posterior probability distributions from the right causal models

Code 
param_grid <- read.csv(file = file.path(path, "direct_and_total_effect_model_data_frame.csv"))
param_grid <- param_grid[param_grid$effect == "direct", ]

param_grid$file <- paste0(remove_special_chars(param_grid$model), ".RDS")

rdss_24 <- list.files(file.path(path, "averaged_models"), pattern = "24.RDS", full.names = TRUE)

combined_draws_list <- lapply(rdss_24, function(x) {

    total_draws <- readRDS(x)

    exp <- attributes((attributes(total_draws)$averaged_fit_formulas))$exposure.formula
    exp <- gsub("n_call | resp_rate(rec_time)  ~  ", "", exp, fixed = TRUE)
    exposure <- exp <- gsub(" + ar(p = 2, time = hour_diff, gr = hour)", "", exp,
        fixed = TRUE)
    exp <- paste0("b_", exp)
    total_draws <- total_draws[, colnames(total_draws) == exp, drop = FALSE]
    names(total_draws) <- exp

    direct_fit_file <- param_grid$file[param_grid$exposure == exposure]

    direct_fit_file <- direct_fit_file[!duplicated(direct_fit_file)]

    if (length(direct_fit_file) > 1)
        direct_fit_file <- grep("24", direct_fit_file, value = TRUE)

    direct_fit <- readRDS(file = file.path(path, "averaged_models", direct_fit_file))
    direct_draws <- posterior::merge_chains(as_draws(direct_fit, variable = exp))
    direct_draws <- as.data.frame(thin_draws(direct_draws, thin = length(direct_draws[[1]][[1]])/(nrow(total_draws)))[[1]])

    direct_draws$effect <- "direct"
    total_draws$effect <- "total"

    draws <- rbind(direct_draws, total_draws)

    return(draws)
})

combined_draws <- do.call(cbind, combined_draws_list)
combined_draws <- combined_draws[, c(which(sapply(combined_draws, is.numeric)), ncol(combined_draws))]

combined_draws[, -ncol(combined_draws)] <- to_change_percentage(combined_draws[,
    -ncol(combined_draws)])


# combined_draws <- as.data.frame(combined_draws)
combined_draws$effect <- ifelse(combined_draws$effect == "direct", "Direct", "Total")


saveRDS(combined_draws, file.path(path, "combined_draws_for_total_and_direct_effects_24h_previous_rain.RDS"))
Code 
combined_draws <- readRDS(file.path(path, "combined_draws_for_total_and_direct_effects_24h_previous_rain.RDS"))

fill_colors <- viridis::mako(10)[c(8, 4)]

gg_dists <- draw_extended_summary(draws = combined_draws, highlight = TRUE, remove.intercepts = TRUE,
    fill = adjustcolor(fill_colors, alpha.f = 0.4), by = "effect", gsub.pattern = c("b_sc_HR",
        "b_sc_rain$", "b_sc_rain_24", "b_sc_temp", "b_sc_moonlight"), gsub.replacement = c("Relative\nhumidity",
        "Current\nrain", "Prior\nrain", "Temperature", "Moonlight"), ylab = "Variable",
    xlab = "Effect size (% of change)")
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Direct-Current rain 3.652 -0.002 7.942 1.009 999.221 949.102
Total-Current rain 0.000 -3.373 3.596 1.001 942.790 942.018
Direct-Moonlight -12.018 -17.074 -6.217 1.001 972.135 944.341
Total-Moonlight -9.600 -15.624 -3.860 1.04 40.665 750.448
Direct-Prior rain 10.060 5.793 14.735 1 916.714 983.283
Total-Prior rain 9.189 5.189 13.581 1.01 883.631 733.445
Direct-Relative humidity 37.000 27.013 47.025 1.001 902.884 848.891
Total-Relative humidity 36.187 27.104 46.613 1.006 942.064 903.519
Direct-Temperature 77.591 65.145 92.509 1.002 955.952 875.689
Total-Temperature 36.112 29.497 43.428 1.043 44.190 893.839

Posterior distribution of direct (green) and total (purple) effect sizes of environmental factors on the calling activity of A. lemur. Posterior values were transformed into percentage change to facilitate interpretation. Dots and error bars show the median and 95% uncertainty intervals of the distributions. Solid color distributions correspond to effect sizes in which uncertainty intervals do not include zero. “Prior rain” accounts for the 24 hour period before sampling calling activity.
Code 
gg_dists + ggplot2::scale_fill_manual(values = fill_colors) + ggplot2::theme(axis.ticks.length = ggplot2::unit(0,
    "pt"), plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"), legend.position = "inside",
    legend.position.inside = c(0.7, 0.7))

Posterior distribution of direct (green) and total (purple) effect sizes of environmental factors on the calling activity of A. lemur. Posterior values were transformed into percentage change to facilitate interpretation. Dots and error bars show the median and 95% uncertainty intervals of the distributions. Solid color distributions correspond to effect sizes in which uncertainty intervals do not include zero. “Prior rain” accounts for the 24 hour period before sampling calling activity.
Code 
param_grid <- read.csv(file = file.path(path, "direct_and_total_effect_model_data_frame.csv"))
param_grid <- param_grid[param_grid$effect == "direct", ]

param_grid$file <- paste0(remove_special_chars(param_grid$model), ".RDS")

rdss_48 <- list.files(file.path(path, "averaged_models"), pattern = "48.RDS", full.names = TRUE)

combined_draws_list <- lapply(rdss_48, function(x) {

    total_draws <- readRDS(x)

    exp <- attributes((attributes(total_draws)$averaged_fit_formulas))$exposure.formula
    exp <- gsub("n_call | resp_rate(rec_time)  ~  ", "", exp, fixed = TRUE)
    exposure <- exp <- gsub(" + ar(p = 2, time = hour_diff, gr = hour)", "", exp,
        fixed = TRUE)
    exp <- paste0("b_", exp)
    total_draws <- total_draws[, colnames(total_draws) == exp, drop = FALSE]
    names(total_draws) <- exp

    direct_fit_file <- param_grid$file[param_grid$exposure == exposure]

    direct_fit_file <- direct_fit_file[!duplicated(direct_fit_file)]

    if (length(direct_fit_file) > 1)
        direct_fit_file <- grep("48", direct_fit_file, value = TRUE)

    direct_fit <- readRDS(file = file.path(path, "averaged_models", direct_fit_file))
    direct_draws <- posterior::merge_chains(as_draws(direct_fit, variable = exp))
    direct_draws <- as.data.frame(thin_draws(direct_draws, thin = length(direct_draws[[1]][[1]])/(nrow(total_draws)))[[1]])

    direct_draws$effect <- "direct"
    total_draws$effect <- "total"

    draws <- rbind(direct_draws, total_draws)

    return(draws)
})

combined_draws <- do.call(cbind, combined_draws_list)
combined_draws <- combined_draws[, c(which(sapply(combined_draws, is.numeric)), ncol(combined_draws))]

combined_draws[, -ncol(combined_draws)] <- to_change_percentage(combined_draws[,
    -ncol(combined_draws)])


# combined_draws <- as.data.frame(combined_draws)
combined_draws$effect <- ifelse(combined_draws$effect == "direct", "Direct", "Total")

saveRDS(combined_draws, file.path(path, "combined_draws_for_total_and_direct_effects_48h_previous_rain.RDS"))
Code 
combined_draws <- readRDS(file.path(path, "combined_draws_for_total_and_direct_effects_48h_previous_rain.RDS"))


gg_dists <- draw_extended_summary(draws = combined_draws, highlight = TRUE, remove.intercepts = TRUE,
    fill = adjustcolor(fill_colors, alpha.f = 0.4), by = "effect", gsub.pattern = c("b_sc_HR",
        "b_sc_rain$", "b_sc_rain_48", "b_sc_temp", "b_sc_moonlight"), gsub.replacement = c("Relative\nhumidity",
        "Current\nrain", "Prior\nrain", "Temperature", "Moonlight"), ylab = "Variable",
    xlab = "Effect size (% of change)")
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Direct-Current rain 2.940 -0.477 6.832 1 933.292 847.036
Total-Current rain -0.544 -4.093 3.071 1 1024.455 819.064
Direct-Moonlight -12.018 -17.074 -6.217 1.001 972.135 944.341
Total-Moonlight -9.925 -15.276 -3.765 1.02 111.506 659.054
Direct-Prior rain 0.623 -2.979 4.124 1.003 857.680 1009.747
Total-Prior rain 0.540 -3.268 4.427 1.014 711.511 873.701
Direct-Relative humidity 36.477 26.300 46.337 0.999 1043.482 878.337
Total-Relative humidity 35.405 25.805 44.967 1.006 801.527 931.946
Direct-Temperature 76.199 63.368 89.305 1.004 1009.050 835.503
Total-Temperature 35.603 29.387 42.312 1.005 961.292 983.283

Code 
gg_dists + ggplot2::scale_fill_manual(values = fill_colors) + ggplot2::theme(axis.ticks.length = ggplot2::unit(0,
    "pt"), plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"), legend.position = "inside",
    legend.position.inside = c(0.7, 0.7))

Takeaways

  • Variation in call activity strongly linked to environmental variation

Session information

R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Costa_Rica
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] pbapply_1.7-2      tidybayes_3.0.6    ggdag_0.2.12       dagitty_0.3-4     
 [5] ohun_1.0.2         warbleR_1.1.30     NatureSounds_1.0.4 seewave_2.2.3     
 [9] tuneR_1.4.6        brmsish_1.0.0      lunar_0.2-1        ggplot2_3.5.1     
[13] ggdist_3.3.2       knitr_1.47         kableExtra_1.4.0   HDInterval_0.2.4  
[17] readxl_1.4.3       posterior_1.5.0    cowplot_1.1.3      brms_2.21.0       
[21] Rcpp_1.0.12        viridis_0.6.5      viridisLite_0.4.2  remotes_2.5.0     
[25] sketchy_1.0.3     

loaded via a namespace (and not attached):
  [1] tensorA_0.36.2.1     rstudioapi_0.16.0    jsonlite_1.8.8      
  [4] magrittr_2.0.3       TH.data_1.1-2        estimability_1.4.1  
  [7] farver_2.1.2         rmarkdown_2.27       vctrs_0.6.5         
 [10] memoise_2.0.1        RCurl_1.98-1.14      htmltools_0.5.8.1   
 [13] distributional_0.4.0 curl_5.2.0           signal_1.8-0        
 [16] cellranger_1.1.0     StanHeaders_2.32.9   KernSmooth_2.23-20  
 [19] htmlwidgets_1.6.4    plyr_1.8.9           sandwich_3.1-0      
 [22] testthat_3.2.1.1     cachem_1.1.0         emmeans_1.9.0       
 [25] zoo_1.8-12           igraph_2.0.3         lifecycle_1.0.4     
 [28] pkgconfig_2.0.3      Matrix_1.6-5         R6_2.5.1            
 [31] fastmap_1.2.0        digest_0.6.36        colorspace_2.1-0    
 [34] labeling_0.4.3       fansi_1.0.6          polyclip_1.10-6     
 [37] abind_1.4-5          compiler_4.3.2       proxy_0.4-27        
 [40] withr_3.0.0          backports_1.5.0      inline_0.3.19       
 [43] DBI_1.2.2            QuickJSR_1.2.2       pkgbuild_1.4.4      
 [46] highr_0.11           ggforce_0.4.2        MASS_7.3-55         
 [49] rjson_0.2.21         classInt_0.4-10      loo_2.7.0           
 [52] tools_4.3.2          units_0.8-5          ape_5.8             
 [55] fftw_1.0-8           glue_1.7.0           nlme_3.1-155        
 [58] grid_4.3.2           sf_1.0-15            checkmate_2.3.1     
 [61] reshape2_1.4.4       generics_0.1.3       diffobj_0.3.5       
 [64] gtable_0.3.5         class_7.3-20         tidyr_1.3.1         
 [67] tidygraph_1.3.1      xml2_1.3.6           utf8_1.2.4          
 [70] ggrepel_0.9.5        pillar_1.9.0         stringr_1.5.1       
 [73] splines_4.3.2        tweenr_2.0.2         dplyr_1.1.4         
 [76] lattice_0.20-45      survival_3.2-13      tidyselect_1.2.1    
 [79] arrayhelpers_1.1-0   gridExtra_2.3        V8_4.4.2            
 [82] svglite_2.1.3        stats4_4.3.2         xfun_0.45           
 [85] graphlayouts_1.1.1   bridgesampling_1.1-2 brio_1.1.4          
 [88] matrixStats_1.3.0    rstan_2.32.6         stringi_1.8.4       
 [91] yaml_2.3.8           boot_1.3-28          xaringanExtra_0.7.0 
 [94] evaluate_0.24.0      codetools_0.2-18     dtw_1.23-1          
 [97] ggraph_2.2.1         tibble_3.2.1         cli_3.6.3           
[100] RcppParallel_5.1.7   xtable_1.8-4         systemfonts_1.1.0   
[103] munsell_0.5.1        coda_0.19-4.1        svUnit_1.0.6        
[106] parallel_4.3.2       rstantools_2.4.0     bayesplot_1.11.1    
[109] Brobdingnag_1.2-9    bitops_1.0-7         mvtnorm_1.2-5       
[112] scales_1.3.0         e1071_1.7-14         packrat_0.9.2       
[115] purrr_1.0.2          crayon_1.5.2         rlang_1.1.4         
[118] multcomp_1.4-25      formatR_1.14